COVID19 UK

Prepare data

Read and format data

Prevalence

df_uk_prev <- read_csv('UK_timeseries_prep_2005.csv')

Parsed with column specification:
cols(
  ut_area = col_character(),
  date = col_character(),
  cumcase = col_double(),
  poptotal = col_double(),
  rate = col_double()
)

df_uk_prev <- df_uk_prev %>% 
  select(ut_area, date, rate) %>% 
  rename(rate_day = rate) %>%
  mutate(date = as.Date(date, "%d%b%Y"))

df_uk_prev

Personality

df_uk_pers <- read_csv('timeseries_uk_utla_march9_april_09.csv')

Parsed with column specification:
cols(
  ut_area = col_character(),
  time = col_double(),
  areaname = col_character(),
  open = col_double(),
  extra = col_double(),
  agree = col_double(),
  neuro = col_double(),
  sci = col_double(),
  frequ = col_double(),
  ut_name = col_character(),
  poptotal = col_double(),
  rate_day = col_double()
)

df_uk_pers <- df_uk_pers %>% 
  select(ut_area, open, agree, neuro, sci, extra) %>% 
  dplyr::rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct()

df_uk_pers

NA

Social distancing

df_uk_socdist <- read_csv('UK_socdist_fb_nuts3.csv')

Parsed with column specification:
cols(
  nuts3 = col_character(),
  date = col_date(format = ""),
  all_day_bing_tiles_visited_relat = col_double(),
  all_day_ratio_single_tile_users = col_double(),
  open = col_double(),
  extra = col_double(),
  agree = col_double(),
  neuro = col_double(),
  sci = col_double(),
  frequ = col_double(),
  nuts3_name = col_character(),
  runday = col_double()
)

df_uk_socdist$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-03-01" "2020-03-08" "2020-03-16" "2020-03-16" "2020-03-24" "2020-03-31" 

df_uk_socdist <- df_uk_socdist %>% select(-runday, -frequ) %>%
  dplyr::rename(pers_o = open, 
                pers_c = sci,
                pers_e = extra,
                pers_a = agree,
                pers_n = neuro) %>% 
  select(-nuts3_name) %>% 
  dplyr::rename(socdist_tiles = all_day_bing_tiles_visited_relat,
                socdist_single_tile = all_day_ratio_single_tile_users) %>%
  drop_na()

df_uk_socdist

Controls

df_uk_ctrl_nuts <- read_csv("controls_UK_nuts3.csv")

Parsed with column specification:
cols(
  nuts3 = col_character(),
  nuts3_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_nuts <- df_uk_ctrl_nuts %>% select(-nuts3_name)
df_uk_ctrl_nuts

df_uk_ctrl_ut <- read_csv("controls_UK_ut.csv")

Parsed with column specification:
cols(
  ut_area = col_character(),
  ut_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_ut <- df_uk_ctrl_ut %>% select(-ut_name)
df_uk_ctrl_ut

NA
NA

Merge prevalence data

df_uk_prev <- df_uk_prev %>% 
  plyr::join(df_uk_pers, by='ut_area') %>% 
  plyr::join(df_uk_ctrl_ut, by='ut_area')

# create sequence of dates
date_sequence <- seq.Date(min(df_uk_prev$date),
                          max(df_uk_prev$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge day index with gps data
df_uk_prev = df_uk_prev %>% 
  merge(df_dates, by='date') %>% 
  arrange(ut_area) %>%
  as_tibble()

df_uk_prev

Merge social distancing data

nuts_ut_key <- read_csv('nuts3_ut.csv')

Parsed with column specification:
cols(
  nuts3 = col_character(),
  ut_area = col_character()
)

df_uk_socdist <- df_uk_socdist %>% plyr::join(df_uk_ctrl_nuts, by='nuts3')

df_uk_socdist <- nuts_ut_key %>% 
  inner_join(df_uk_socdist, by = c('nuts3')) %>%
  inner_join(select(df_uk_prev, ut_area, date, rate_day), by = c('ut_area', 'date')) %>%
  select(-ut_area)

# create sequence of dates
date_sequence <- seq.Date(min(df_uk_socdist$date),
                          max(df_uk_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge day index with gps data
df_uk_socdist = df_uk_socdist %>% 
  merge(df_dates, by='date') %>% 
  arrange(nuts3) %>%
  as_tibble()


df_uk_socdist

NA

Identify London areas

nuts_london_inner <- c('UKI31','UKI32','UKI33','UKI34','UKI41',
                      'UKI42','UKI43','UKI44','UKI45')

nuts_london_outer <- c('UKI51','UKI52','UKI53','UKI54','UKI61',
                      'UKI62','UKI63','UKI71','UKI72','UKI73',
                      'UKI74','UKI75')

ut_london_inner <- c('E09000007','E09000001','E09000033','E09000013',
                    'E09000020','E09000032','E09000025','E09000012',
                    'E09000030','E09000014','E09000019','E09000023',
                    'E09000028','E09000022')

ut_london_outer <- c('E09000011','E09000004','E09000016','E09000002',
                    'E09000031','E09000026','E09000010','E09000006',
                    'E09000008','E09000029','E09000021','E09000024',
                    'E09000003','E09000005','E09000009','E09000017',
                    'E09000015','E09000018','E09000027')

df_uk_prev = df_uk_prev %>% 
  mutate(london = ifelse(ut_area %in% ut_london_inner, 'london_inner', 
                       ifelse(ut_area %in% ut_london_outer, 'london_outer',
                              'country'))) %>%
  mutate(london = as.factor(london))

df_uk_socdist = df_uk_socdist %>% 
  mutate(london = ifelse(nuts3 %in% nuts_london_inner, 'london_inner', 
                       ifelse(nuts3 %in% nuts_london_outer, 'london_outer',
                              'country'))) %>%
  mutate(london = as.factor(london))

Check timeframes

df_uk_prev$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-01-30" "2020-02-26" "2020-03-24" "2020-03-24" "2020-04-21" "2020-05-18" 

df_uk_socdist$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-03-01" "2020-03-08" "2020-03-16" "2020-03-16" "2020-03-24" "2020-03-31" 

Control for weekend effect in social distancing

df_uk_loess <- df_uk_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!weekday %in% c('6','7')) %>% 
  split(.$nuts3) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_uk_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'nuts3', value = 'loess') %>% 
  group_by(nuts3) %>% 
  mutate(time = row_number())

df_uk_loess

df_uk_socdist <- df_uk_socdist %>% merge(df_uk_loess, by=c('nuts3', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7'), loess,
                                            socdist_single_tile)) %>%
  arrange(nuts3, time) %>% 
  select(-weekday)


df_uk_socdist <- df_uk_socdist %>% drop_na() %>% mutate(time = time-1)

Explore data

Plot prevalence over time

df_uk_prev %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_prev %>% mutate(prev_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Explore differences between london and the rest

df_uk_prev %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  facet_wrap(~london) +
  ggtitle("Overall prevalence over time")

NA
NA
NA

Plot social distancing over time

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  facet_wrap(~london) +
  ggtitle("Overall social distancing over time")

df_uk_socdist <- df_uk_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>% 
  select(-loess, -socdist_single_tile_clean)

Correlations

df_uk_prev %>% group_by(ut_area) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-ut_area, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3) %>%
  as.data.frame()

df_uk_socdist %>% group_by(nuts3) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-nuts3, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3) %>% 
  as.data.frame()

NA

Rescale Data

lvl2_scaled_ut <- df_uk_prev %>% 
  dplyr::select(-time, -date, -rate_day, -london) %>% 
  distinct() %>% 
  mutate_at(vars(-ut_area), scale)

lvl1_scaled_ut <- df_uk_prev %>% select(ut_area, time, rate_day)

df_uk_prev_scaled <- plyr::join(lvl1_scaled_ut, lvl2_scaled_ut, by = 'ut_area')

df_uk_prev_scaled

lvl2_scaled_nuts <- df_uk_socdist %>% 
  dplyr::select(-time, -date, -london, 
                -socdist_tiles, -socdist_single_tile, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-nuts3), scale)

lvl1_scaled_nuts <- df_uk_socdist %>% select(nuts3, time, socdist_single_tile, rate_day) %>% 
  mutate_at(vars(-nuts3, -time, - rate_day), scale)

df_uk_socdist_scaled <- plyr::join(lvl1_scaled_nuts, lvl2_scaled_nuts, by = 'nuts3')

df_uk_socdist_scaled

NA

Predict Prevalence

Extract first day of covid outbreak

# get onset day
df_uk_onset_prev <- df_uk_prev_scaled %>% 
  group_by(ut_area) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time))
  
# merge with county data
df_uk_onset_prev <- df_uk_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  left_join(df_uk_onset_prev, by = 'ut_area')

# handle censored data
df_uk_onset_prev <- df_uk_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_uk_prev$date)))+1))

df_uk_onset_prev

NA

Extract slopes

# cut time series before onset
df_uk_prev_scaled <- df_uk_prev_scaled %>% 
  group_by(ut_area) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_uk_prev_scaled <- df_uk_prev_scaled %>%
  group_by(ut_area) %>%
  filter(n() == 30) %>%
  ungroup()

# extract slope prevalence
df_uk_slope_prev <- df_uk_prev_scaled %>% split(.$ut_area) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('ut_area') %>% 
  rename(slope_prev = '.')

# merge with county data
df_uk_slope_prev <- df_uk_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_uk_slope_prev, by = 'ut_area') %>%
  drop_na()

Explore distributions

df_uk_onset_prev %>% ggplot(aes(onset_prev)) + geom_histogram()

df_uk_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram()

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_uk_onset_prev)
cox_onset_prev %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_uk_onset_prev)

  n= 149, number of events= 149 

            coef exp(coef)  se(coef)      z Pr(>|z|)  
pers_o  0.322264  1.380249  0.166864  1.931   0.0534 .
pers_c  0.036475  1.037149  0.135119  0.270   0.7872  
pers_e -0.013895  0.986201  0.143396 -0.097   0.9228  
pers_a  0.005323  1.005338  0.106927  0.050   0.9603  
pers_n -0.092400  0.911740  0.106197 -0.870   0.3843  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.3802     0.7245    0.9952     1.914
pers_c    1.0371     0.9642    0.7958     1.352
pers_e    0.9862     1.0140    0.7446     1.306
pers_a    1.0053     0.9947    0.8153     1.240
pers_n    0.9117     1.0968    0.7404     1.123

Concordance= 0.605  (se = 0.027 )
Likelihood ratio test= 12.64  on 5 df,   p=0.03
Wald test            = 13.47  on 5 df,   p=0.02
Score (logrank) test = 13.78  on 5 df,   p=0.02

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                             data = df_uk_onset_prev)
cox_onset_prev_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative, 
    data = df_uk_onset_prev)

  n= 144, number of events= 144 
   (5 observations deleted due to missingness)

                  coef exp(coef) se(coef)      z Pr(>|z|)  
pers_o         0.13431   1.14375  0.26537  0.506   0.6128  
pers_c         0.27571   1.31746  0.21410  1.288   0.1978  
pers_e        -0.28546   0.75167  0.16166 -1.766   0.0774 .
pers_a         0.15618   1.16904  0.12936  1.207   0.2273  
pers_n        -0.04289   0.95801  0.13821 -0.310   0.7563  
airport_dist   0.08819   1.09219  0.10241  0.861   0.3892  
males         -0.28751   0.75013  0.11525 -2.495   0.0126 *
popdens        0.31764   1.37388  0.22294  1.425   0.1542  
manufacturing -0.01434   0.98576  0.15053 -0.095   0.9241  
tourism        0.20101   1.22263  0.12897  1.559   0.1191  
health        -0.13046   0.87769  0.09930 -1.314   0.1889  
academic       0.41323   1.51170  0.25879  1.597   0.1103  
medinc        -0.10312   0.90201  0.15281 -0.675   0.4998  
medage        -0.54104   0.58214  0.21413 -2.527   0.0115 *
conservative  -0.05175   0.94957  0.28447 -0.182   0.8557  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.1438     0.8743    0.6799    1.9241
pers_c           1.3175     0.7590    0.8660    2.0044
pers_e           0.7517     1.3304    0.5475    1.0319
pers_a           1.1690     0.8554    0.9072    1.5064
pers_n           0.9580     1.0438    0.7307    1.2561
airport_dist     1.0922     0.9156    0.8936    1.3350
males            0.7501     1.3331    0.5985    0.9402
popdens          1.3739     0.7279    0.8875    2.1267
manufacturing    0.9858     1.0144    0.7339    1.3240
tourism          1.2226     0.8179    0.9496    1.5742
health           0.8777     1.1394    0.7225    1.0663
academic         1.5117     0.6615    0.9103    2.5104
medinc           0.9020     1.1086    0.6686    1.2170
medage           0.5821     1.7178    0.3826    0.8857
conservative     0.9496     1.0531    0.5437    1.6583

Concordance= 0.643  (se = 0.027 )
Likelihood ratio test= 37.26  on 15 df,   p=0.001
Wald test            = 39.13  on 15 df,   p=6e-04
Score (logrank) test = 40.98  on 15 df,   p=3e-04

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_uk_slope_prev)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.032324 -0.011515 -0.003789  0.012728  0.048348 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0239193  0.0014399  16.612  < 2e-16 ***
pers_o      -0.0060442  0.0028372  -2.130  0.03492 *  
pers_c      -0.0074351  0.0023785  -3.126  0.00216 ** 
pers_e       0.0033599  0.0026519   1.267  0.20731    
pers_a      -0.0014460  0.0020483  -0.706  0.48139    
pers_n      -0.0001847  0.0021315  -0.087  0.93109    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01727 on 138 degrees of freedom
Multiple R-squared:  0.121, Adjusted R-squared:  0.0892 
F-statistic: 3.801 on 5 and 138 DF,  p-value: 0.002952

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                         data = df_uk_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative, 
    data = df_uk_slope_prev)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.028137 -0.010922 -0.002334  0.009211  0.045494 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    2.388e-02  1.380e-03  17.303   <2e-16 ***
pers_o        -3.760e-03  4.133e-03  -0.910   0.3648    
pers_c        -7.549e-03  3.233e-03  -2.335   0.0211 *  
pers_e         1.386e-03  2.733e-03   0.507   0.6128    
pers_a         4.417e-04  2.399e-03   0.184   0.8542    
pers_n        -1.024e-03  2.440e-03  -0.419   0.6756    
airport_dist  -4.622e-03  1.818e-03  -2.542   0.0122 *  
males         -1.683e-03  2.049e-03  -0.821   0.4131    
popdens        6.486e-03  3.510e-03   1.848   0.0669 .  
manufacturing  8.049e-05  2.368e-03   0.034   0.9729    
tourism       -8.124e-04  2.007e-03  -0.405   0.6863    
health         1.847e-03  1.739e-03   1.062   0.2904    
academic      -8.043e-04  4.209e-03  -0.191   0.8488    
medinc         3.966e-03  2.549e-03   1.556   0.1221    
medage         4.778e-03  3.375e-03   1.416   0.1593    
conservative   5.636e-03  4.576e-03   1.231   0.2204    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01653 on 128 degrees of freedom
Multiple R-squared:  0.253, Adjusted R-squared:  0.1655 
F-statistic:  2.89 on 15 and 128 DF,  p-value: 0.0006055

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                           data = df_uk_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
nuts_complete <- df_uk_socdist_scaled %>% 
  group_by(nuts3) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$nuts3

# run changepoint analysis
df_uk_socdist_cpt_results <- df_uk_socdist_scaled %>% 
  select(nuts3, socdist_single_tile) %>%
  filter(nuts3 %in% nuts_complete) %>% 
  split(.$nuts3) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_uk_socdist_cpt_day <- df_uk_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('nuts3')

# calculate mean differences
df_uk_socdist_cpt_mean_diff <- df_uk_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('nuts3')

# calculate varaince differences
df_uk_socdist_cpt_var_diff <- df_uk_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$variance) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(var_diff_socdist = '.') %>%
  rownames_to_column('nuts3')

# merge with county data
df_uk_cpt_socdist <- df_uk_socdist_scaled %>% 
  select(-time, -rate_day, -socdist_single_tile) %>%
  distinct() %>% 
  left_join(df_uk_socdist_cpt_day, by='nuts3') %>%
  left_join(df_uk_socdist_cpt_mean_diff, by='nuts3') %>%
  left_join(df_uk_socdist_cpt_var_diff, by='nuts3') %>%
  left_join(nuts_ut_key, by='nuts3') %>% 
  left_join(select(df_uk_onset_prev, ut_area, onset_prev), by='ut_area') %>%
  left_join(select(df_uk_slope_prev, ut_area, slope_prev), by='ut_area') %>%
  select(-ut_area)

# handle censored data
df_uk_cpt_socdist <- df_uk_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), 
                                  as.numeric(diff(range(df_us$date))), 
                                  cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= 60, 0, 1))

df_uk_cpt_socdist$cpt_day_socdist %>% hist()

df_uk_cpt_socdist$mean_diff_socdist %>% hist()

df_uk_cpt_socdist$var_diff_socdist %>% hist()

for(i in head(df_uk_socdist_cpt_results, 5)){
  plot(i)
}

NA

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_uk_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_uk_cpt_socdist)

  n= 131, number of events= 131 

           coef exp(coef) se(coef)      z Pr(>|z|)
pers_o  0.02292   1.02319  0.14906  0.154    0.878
pers_c  0.01178   1.01185  0.14322  0.082    0.934
pers_e  0.17493   1.19117  0.14723  1.188    0.235
pers_a -0.19324   0.82428  0.11907 -1.623    0.105
pers_n -0.01502   0.98509  0.12389 -0.121    0.903

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.0232     0.9773    0.7640     1.370
pers_c    1.0118     0.9883    0.7642     1.340
pers_e    1.1912     0.8395    0.8926     1.590
pers_a    0.8243     1.2132    0.6527     1.041
pers_n    0.9851     1.0151    0.7727     1.256

Concordance= 0.672  (se = 0.046 )
Likelihood ratio test= 11.4  on 5 df,   p=0.04
Wald test            = 13.38  on 5 df,   p=0.02
Score (logrank) test = 13.55  on 5 df,   p=0.02

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                  data = df_uk_cpt_socdist)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_dist + males + popdens + 
    manufacturing + tourism + health + academic + medinc + medage + 
    conservative, data = df_uk_cpt_socdist)

  n= 131, number of events= 131 

                   coef exp(coef)  se(coef)      z Pr(>|z|)  
pers_o        -0.100979  0.903952  0.243690 -0.414   0.6786  
pers_c        -0.048635  0.952528  0.198550 -0.245   0.8065  
pers_e         0.080996  1.084366  0.163732  0.495   0.6208  
pers_a        -0.006914  0.993110  0.162410 -0.043   0.9660  
pers_n         0.204498  1.226909  0.151016  1.354   0.1757  
airport_dist   0.251102  1.285441  0.126736  1.981   0.0476 *
males         -0.148267  0.862201  0.142494 -1.041   0.2981  
popdens        0.158406  1.171641  0.233680  0.678   0.4979  
manufacturing -0.083042  0.920313  0.130937 -0.634   0.5259  
tourism       -0.254837  0.775043  0.142800 -1.785   0.0743 .
health        -0.266804  0.765823  0.122195 -2.183   0.0290 *
academic       0.153738  1.166185  0.247465  0.621   0.5344  
medinc         0.166816  1.181537  0.177460  0.940   0.3472  
medage        -0.008602  0.991435  0.216020 -0.040   0.9682  
conservative  -0.180652  0.834726  0.268518 -0.673   0.5011  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           0.9040     1.1063    0.5607    1.4574
pers_c           0.9525     1.0498    0.6455    1.4057
pers_e           1.0844     0.9222    0.7867    1.4947
pers_a           0.9931     1.0069    0.7224    1.3653
pers_n           1.2269     0.8151    0.9126    1.6495
airport_dist     1.2854     0.7779    1.0027    1.6479
males            0.8622     1.1598    0.6521    1.1400
popdens          1.1716     0.8535    0.7411    1.8523
manufacturing    0.9203     1.0866    0.7120    1.1896
tourism          0.7750     1.2903    0.5858    1.0254
health           0.7658     1.3058    0.6027    0.9731
academic         1.1662     0.8575    0.7180    1.8941
medinc           1.1815     0.8464    0.8344    1.6730
medage           0.9914     1.0086    0.6492    1.5141
conservative     0.8347     1.1980    0.4932    1.4129

Concordance= 0.756  (se = 0.043 )
Likelihood ratio test= 30.64  on 15 df,   p=0.01
Wald test            = 31.24  on 15 df,   p=0.008
Score (logrank) test = 32.93  on 15 df,   p=0.005

# predict hazard from personality with controls
cox_cpt_socdist_ctrl2 <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                  data = df_uk_cpt_socdist)
cox_cpt_socdist_ctrl2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + airport_dist + males + popdens + 
    manufacturing + tourism + health + academic + medinc + medage + 
    conservative + onset_prev + slope_prev, data = df_uk_cpt_socdist)

  n= 130, number of events= 130 
   (1 observation deleted due to missingness)

                    coef  exp(coef)   se(coef)      z Pr(>|z|)  
pers_o        -1.002e-01  9.046e-01  2.493e-01 -0.402   0.6877  
pers_c        -1.163e-01  8.902e-01  2.070e-01 -0.562   0.5742  
pers_e         8.183e-02  1.085e+00  1.631e-01  0.502   0.6158  
pers_a         5.279e-02  1.054e+00  1.731e-01  0.305   0.7604  
pers_n         2.313e-01  1.260e+00  1.571e-01  1.472   0.1411  
airport_dist   2.164e-01  1.242e+00  1.291e-01  1.676   0.0936 .
males         -2.161e-01  8.057e-01  1.560e-01 -1.385   0.1661  
popdens        3.186e-01  1.375e+00  2.541e-01  1.254   0.2098  
manufacturing -6.886e-02  9.335e-01  1.314e-01 -0.524   0.6003  
tourism       -2.531e-01  7.764e-01  1.447e-01 -1.749   0.0803 .
health        -2.746e-01  7.599e-01  1.267e-01 -2.167   0.0302 *
academic       7.433e-02  1.077e+00  2.554e-01  0.291   0.7711  
medinc         2.756e-01  1.317e+00  1.876e-01  1.469   0.1419  
medage         4.406e-02  1.045e+00  2.175e-01  0.203   0.8395  
conservative  -1.681e-01  8.453e-01  2.693e-01 -0.624   0.5324  
onset_prev    -1.956e-03  9.980e-01  1.470e-02 -0.133   0.8941  
slope_prev    -1.053e+01  2.674e-05  9.175e+00 -1.148   0.2511  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o        9.046e-01  1.105e+00 5.550e-01    1.4745
pers_c        8.902e-01  1.123e+00 5.933e-01    1.3356
pers_e        1.085e+00  9.214e-01 7.884e-01    1.4940
pers_a        1.054e+00  9.486e-01 7.509e-01    1.4801
pers_n        1.260e+00  7.935e-01 9.261e-01    1.7148
airport_dist  1.242e+00  8.054e-01 9.641e-01    1.5990
males         8.057e-01  1.241e+00 5.934e-01    1.0939
popdens       1.375e+00  7.271e-01 8.358e-01    2.2628
manufacturing 9.335e-01  1.071e+00 7.215e-01    1.2077
tourism       7.764e-01  1.288e+00 5.846e-01    1.0310
health        7.599e-01  1.316e+00 5.927e-01    0.9741
academic      1.077e+00  9.284e-01 6.529e-01    1.7771
medinc        1.317e+00  7.591e-01 9.120e-01    1.9027
medage        1.045e+00  9.569e-01 6.823e-01    1.6007
conservative  8.453e-01  1.183e+00 4.986e-01    1.4328
onset_prev    9.980e-01  1.002e+00 9.697e-01    1.0272
slope_prev    2.674e-05  3.740e+04 4.145e-13 1724.8956

Concordance= 0.745  (se = 0.044 )
Likelihood ratio test= 32.6  on 17 df,   p=0.01
Wald test            = 32.53  on 17 df,   p=0.01
Score (logrank) test = 33.98  on 17 df,   p=0.008

Linear models predicting mean differences

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_uk_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_uk_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.236161 -0.073887 -0.001491  0.063529  0.288673 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.617663   0.009637 167.864  < 2e-16 ***
pers_o      -0.004645   0.017588  -0.264 0.792145    
pers_c      -0.039209   0.016089  -2.437 0.016220 *  
pers_e       0.059962   0.017020   3.523 0.000596 ***
pers_a      -0.020685   0.014162  -1.461 0.146648    
pers_n      -0.004006   0.014038  -0.285 0.775854    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1103 on 125 degrees of freedom
Multiple R-squared:  0.4319,    Adjusted R-squared:  0.4092 
F-statistic: 19.01 on 5 and 125 DF,  p-value: 5.03e-14

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative,
                            data = df_uk_cpt_socdist)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative, 
    data = df_uk_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.240455 -0.063142 -0.001388  0.051147  0.206542 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.6176625  0.0075553 214.109  < 2e-16 ***
pers_o        -0.0326846  0.0220339  -1.483   0.1407    
pers_c        -0.0130620  0.0177712  -0.735   0.4638    
pers_e         0.0120006  0.0146614   0.819   0.4148    
pers_a         0.0360106  0.0140066   2.571   0.0114 *  
pers_n         0.0229335  0.0122311   1.875   0.0633 .  
airport_dist  -0.0208327  0.0106569  -1.955   0.0530 .  
males         -0.0256885  0.0125746  -2.043   0.0434 *  
popdens        0.0928652  0.0197043   4.713  6.9e-06 ***
manufacturing -0.0252661  0.0117189  -2.156   0.0332 *  
tourism       -0.0103949  0.0112219  -0.926   0.3562    
health        -0.0005235  0.0099268  -0.053   0.9580    
academic       0.0282827  0.0220672   1.282   0.2025    
medinc         0.0263085  0.0145454   1.809   0.0731 .  
medage        -0.0012929  0.0193857  -0.067   0.9469    
conservative  -0.0049255  0.0247582  -0.199   0.8427    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08647 on 115 degrees of freedom
Multiple R-squared:  0.6787,    Adjusted R-squared:  0.6368 
F-statistic:  16.2 on 15 and 115 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl2 <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                            data = df_uk_cpt_socdist)
lm_meandiff_socdist_ctrl2 %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + airport_dist + males + popdens + manufacturing + 
    tourism + health + academic + medinc + medage + conservative + 
    onset_prev + slope_prev, data = df_uk_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.248333 -0.059558  0.000339  0.056215  0.216795 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.6305244  0.0318952  51.121  < 2e-16 ***
pers_o        -0.0353915  0.0223396  -1.584   0.1160    
pers_c        -0.0090910  0.0179360  -0.507   0.6133    
pers_e         0.0125871  0.0146878   0.857   0.3933    
pers_a         0.0357278  0.0140017   2.552   0.0121 *  
pers_n         0.0218672  0.0122562   1.784   0.0771 .  
airport_dist  -0.0178000  0.0107886  -1.650   0.1018    
males         -0.0179556  0.0133454  -1.345   0.1812    
popdens        0.0859734  0.0206309   4.167 6.09e-05 ***
manufacturing -0.0261575  0.0118566  -2.206   0.0294 *  
tourism       -0.0110823  0.0112664  -0.984   0.3274    
health        -0.0008650  0.0099392  -0.087   0.9308    
academic       0.0283153  0.0222260   1.274   0.2053    
medinc         0.0214989  0.0148569   1.447   0.1507    
medage        -0.0004348  0.0194895  -0.022   0.9822    
conservative  -0.0078513  0.0247978  -0.317   0.7521    
onset_prev    -0.0012086  0.0011796  -1.025   0.3078    
slope_prev     1.2542002  0.6996823   1.793   0.0757 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08639 on 112 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.6857,    Adjusted R-squared:  0.638 
F-statistic: 14.38 on 17 and 112 DF,  p-value: < 2.2e-16

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               airport_dist + males + popdens + manufacturing + 
                               tourism + health + academic + medinc + medage +
                               conservative + onset_prev + slope_prev,
                               data = df_uk_cpt_socdist %>% drop_na(),
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))